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Summary. We compute numerically eigenvalues and eigenfunctions of the quantum 
Hamiltonian that describes the quantum mechanics of a point particle moving freely 
in a particular three-dimensional hyperbolic space of finite volume and investigate 
the distribution of the eigenvalues. 



1 Introduction 

The distribution of the eigenvalues of a quantum Hamiltonian is a central 
subject that is studied in quantum chaos. There are some generally accepted 
conjectures about the nearest-neighbor spacing distributions of the eigenval- 
ues. 

Unless otherwise stated we use the following assumptions: The quantum 
mechanical system is desymmetrized with respect to all its unitary symme- 
tries, and whenever we examine the distribution of the eigenvalues we regard 
them on the scale of the mean level spacings. Moreover, it is generically be- 
lieved that after desymmetrization a generic quantum Hamiltonian possesses 
no degenerate eigenvalues. 

Conjecture 1 (Berry, Tabor '11). If the corresponding classical system is in- 
tegrable, the eigenvalues behave like independent random variables and the 
distribution of the nearest-neighbor spacings is close to the Poisson distribu- 
tion, i.e. there is no level repulsion. 

Conjecture 2 (Bohigas, Giannoni, Schmit '2, 31). If the corresponding classi- 
cal system is chaotic, the eigenvalues are distributed like the eigenvalues of 
hermitian random matrices [1] . The corresponding ensembles depend only on 
the symmetries of the system: 

• For chaotic systems without time-reversal invariance the distribution of 
the eigenvalues should be close to the distribution of the Gaussian Unitary 
Ensemble (GUE) which is characterized by a quadratic level repulsion. 
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• For chaotic systems with time-reversal invariance and integer spin the 
distribution of the eigenvalues should be close to the distribution of the 
Gaussian Orthogonal Ensemble (GOE) which is characterized by a linear 
level repulsion. 

• For chaotic systems with time-reversal invariance and half-integer spin the 
distribution of the eigenvalues should be close to the distribution of the 
Gaussian Symplectic Ensemble (GSE) which is characterized by a quartic 
level repulsion. 

These conjectures are very well confirmed by numerical calculations, but sev- 
eral exceptions are known. Here are two examples: 

Exception 1 The harmonic oscillator is classically integrable, but its spec- 
trum is equidistant. 

Exception 2 The geodesic motion on surfaces with constant negative curva- 
ture provides a prime example for classical chaos. In some cases, however, 
the nearest-neighbor distribution of the eigenvalues of the Laplacian on these 
surfaces appears to be Poissonian. 

"A strange arithmetical structure of chaos" in the case of surfaces of constant 
negative curvature that are generated by arithmetic fundamental groups was 
discovered by Aurich and Steiner [5] , see also Aurich, Bogomolny, and Steiner 
UJ- Deviations from the expected GOE-behaviour in the case of a particu- 
lar arithmetic surface were numerically observed by Bohigas, Giannoni, and 
Schmit |3] and by Aurich and Steiner 7 . Computations coming out in [71 [Hj 
showed, however, that the level statistics on 30 generic (i.e. non-arithmetic) 
surfaces were in nice agreement with the expected random-matrix theory pre- 
diction in accordance with conjecture |5] This has led Bogomolny, Georgeot, 
Giannoni, and Schmit [H], Bolte, Steil, and Steiner ^Oj) and Sarnak to 
introduce the concept of arithmetic quantum chaos. 

Conjecture 3 (Arithmetic Quantum Chaos). On surfaces of constant negative 
curvature that are generated by arithmetic fundamental groups, the distribu- 
tion of the eigenvalues of the quantum Hamiltonian are close to the Poisson 
distribution. Due to level clustering small spacings occur comparably often. 

We compute numerically the eigenvalues and eigenfunctions of the Lapla- 
cian that describes the quantum mechanics of a point particle moving freely 
in the non-integrable three-dimensional hyperbolic space of constant negative 
curvature generated by the Picard group. The Picard group is arithmetic and 
we find that our results are in accordance with the conjecture of arithmetic 
quantum chaos. 

For the definition of an arithmetic group we refer the reader to |12j . 
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2 Preliminaries: The modular group 

For simplicity we first introduce the topology and geometry of the two- 
dimensional surface of constant negative curvature that is generated by the 
modular group P2J . It will then be easy to carry over to the three-dimensional 
space of constant negative curvature that is generated by the Picard group. 
The construction begins with the upper half-plane, 



equipped with the hyperbolic metric of constant negative curvature 



A free particle on the upper half-plane moves along geodesies, which are 
straight lines and semicircles perpendicular to the £-axis, respectively, see 
figure 



Fig. 1. Geodesies in the upper half-plane of constant negative curvature. 

Expressing a point {x, y) S Ji as a complex number z — x+iy, all isometries 
of the hyperbolic metric are given by the group of linear fractional transfor- 
mations, 



H = {(x,y)£R 2 ; y>0}, 



ds 2 = 



dx 2 + dy 2 



y 




x 



az + b 



a, b,c,d £ R, ad — be = 1, 



z i— * 7Z = 



cz + ef 



which is isomorphic to the group of matrices 




up to a common sign of the matrix entries, 



H. Then 



SL(2,M)/{±1} = PSL(2,E). 

In analogy to the concept of a fundamental cell in a regular lattice of 
a crystal we can introduce a fundamental domain of a discrete group r C 
PSL(2,M). 

Definition 1. A fundamental domain of the discrete group r is an open sub- 
set T C Ti with the following conditions: The closure of T meets each orbit 
Fz = {jz; 7 S r} at least once, T meets each orbit Fz at most once, and 
the boundary of T has Lebesgue measure zero. 

If we choose the group r to be the modular group, 

r = PSL(2,Z), 
which is generated by a translation and an inversion, 



1 1 
1 

-1 

1 



z i— > z + 1, 



Z I > — z 



the fundamental domain of standard shape is 

T = {z — x + \y e H; 
see figure |3 The isometric copies of the fundamental domain 7 e r, 



1 1 

2 <x< 2> 1^1 > 1 >' 



-1/2 



1/2 



Fig. 2. The fundamental domain of the modular group. 



tessellate the upper half-plane completely without any overlap or gap, see 
figure 
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Fig. 3. The upper half-plane tessellated with isometric copies of the fundamental 
domain. 

Identifying the fundamental domain T and parts of its boundary with all 
its isometric copies -fj 7 , V7 £ r, defines the topology to be the quotient space 
r\H. The quotient space r\Tt can also be thought of as the fundamental 
domain T with its faces glued according to the elements of the group r, see 
figure 0] 








-1/2 


1/2 


X 



Fig. 4. Identifying the faces of the fundamental domain according to the elements 
of the modular group. 

Any function being defined on the upper half-plane which is invariant 
under linear fractional transformations, 

f(z) = f{yz) V 7 er, 

can be identified with a function living on the quotient space r\H. A function 
on the quotient space is tantamount to a function on the fundamental domain 
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with periodic boundary conditions. Vice versa, any function being defined on 
the quotient space can be identified with an automorphic function, f(z) — 
f{")z), V7 € r, living on the upper half-plane. 

With the hyperbolic metric the quotient space r\H inherits the structure 
of an orbifold. An orbifold locally looks like a manifold, with the exception 
that it is allowed to have elliptic fix-points. 

The orbifold of the modular group has one parabolic and two elliptic fix- 
points, 

a 1 -V3 

2 = 100, z = 1, and z = 2 ~2~' 

The parabolic one fixes a cusp at z = ioo which is invariant under the parabolic 
element 

1 1 
1 

Hence the orbifold of the modular group is non-compact. The volume element 
corresponding to the hyperbolic metric reads 

dxdy 
d/i = — 2-, 

such that the volume of the orbifold r\H is finite, 

vol(AW) = |. 

Scaling the units such that H = 1 and 2m = 1, the stationary Schrodingcr 
equation which describes the quantum mechanics of a point particle moving 
freely in the orbifold r\H becomes 

{A + X)f(z) = 0, 

where the hyperbolic Laplacian is given by 

A = y 2 (— + — 

dx 2 dy 2 

and A is the scaled energy. We can relate the the eigenvalue problem defined on 
the orbifold r\H to the eigenvalue problem defined on the upper-half space, 
with the eigenfunctions being subject to the automorphy condition relative to 
the discrete group _T, 

f(ryz) = f{z) v 7 er. 



In order to avoid solutions that grow exponentially in the cusp, we impose 
the boundary condition 
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f(z) — 0(y K ) for z — > ioo 

where k is some positive constant. 

The solutions of this eigenvalue problem can be identified with Maass 
waveforms (14] . The identification is worthwhile, since much is known about 
Maass waveforms from number theory and harmonic analysis which will sim- 
plify their computation, see e.g. [T31[TCl[r7l[THl[TUlOniOmil^lE5|. 



3 The Picard group 

In the three-dimensional case one considers the upper-half space, 

H = {(x ,x 1 ,y)£m 3 ; y > 0} 
equipped with the hyperbolic metric 

2 _ dxl + dx\ + dy 2 

y 2 

The geodesies of a particle moving freely in the upper half-space are straight 
lines and semicircles perpendicular to the xo-xi-plane, respectively, see figure 

El 




x o 



Fig. 5. Geodesies in the upper half-space of constant negative curvature. 



Expressing any point (xo, Xi, y) £ H as a Hamilton quaternion, z = xq + 
ixi +jy, with the multiplication defined by i 2 = — 1, j 2 = — 1, ij +ji = 0, all 
motions in the upper half-space are given by linear fractional transformations 

z jz — (az + b)(cz + d) , a, b, c, d g C, ad—bc=l. 

The group of these transformations is isomorphic to the group of matrices 



7 = M GSL(2,C) 



H. Then 

up to a common sign of the matrix entries, 

SL(2,C)/{±1} = PSL(2,C). 

The motions provided by the elements of PSL(2, <D) exhaust all orientation 
preserving isometries of the hyperbolic metric on H. 

Remark 1. If one wants to avoid using quaternions, the point (xQ,x\,y) G TL 
can be expressed by (x, y) G <D x R with x — xq + ixi and y > 0. But then 
the linear fractional transformation look somewhat more complicated, 

/ \ / \ / (ax + b) (ex + d) + acy 2 y , 

(x,y) i-> l(x,y) = ( ; — jo , | I — r^J- 

v |cx + d\ z + \cy\ z \cx + d\ z + \cy\" 

In order to keep the notation simple we hence use quaternions. 

We now choose the discrete group r C PSL(2, C) generated by the cosets 
of three elements, 

1 1\ (1 i\ (0 -1 

which yield two translations and one inversion, 

z h- > z + 1, zH>z|i, z i ^ — z . 

This group I 1 is called the Picard group. The three motions generating r, 
together with the coset of the element 

i 
-i 

that is isomorphic to the symmetry 

z = x + jy h-> izi = -x + jy, 
can be used to construct the fundamental domain of standard shape 

T = {z = x + ixi + jy G H; ~ < x < i, < x x < i, \z\ > 1}, 

see figure |H| Identifying the faces of the fundamental domain according to the 
elements of the group r leads to a realization of the quotient space r\Tt, see 
figure H 

With the hyperbolic metric the quotient space r\H inherits the structure 
of an orbifold that has one parabolic and four elliptic fix-points, 
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Fig. 6. The fundamental domain of the Picard group. 




Fig. 7. Identifying the faces of the fundamental domain according to the elements 
of the Picard group. 



1 .1 . 1 

2 +1 2 +J V2' 



.1 . /3 
1 2 +J VI- 



The parabolic fix-point corresponds to a cusp at z = joo that is invariant 
under the parabolic elements 



1 1 
1 



and 



1 i 
1 



The volume element deriving from the hyperbolic metric reads 

dx^dxidy 



dpi = 



y 



such that the volume of the non-compact orbifold r\H is finite 

Ck(2) 



voi(r\w) 



4tt 2 



0.305 
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where 

^ez[i]-{o} 

is the Dedekind zeta function. 

We are interested in the eigenfunctions of the Laplacian, 

2 ,_c^_ a 2 v _9 

^ 9xq dx 2 <9y 2 ' ^ 9y ' 

which determine the quantum mechanics of a particle moving freely in the 
orbifold r\TC. As in the preceding section we identify the solutions with Maass 
waveforms |25) . 

Since the Maass waveforms are automorphic, and therefore periodic in xo 
and x\, it follows that they can be expanded into a Fourier series, 

f(z) = u(y) + Yl a p yK ir (27r\(3\y)e 27rmPx , (1) 

/3ez[i]-{o} 



where 

u(y) 



b y 1+ir + hy 1 -^ if r ^ 0, 
biy + b3ylny if r = 0. 



K- lr (x) is the K-Bcssel function whose order is connected with the eigenvalue 
A by 



A = r 2 + 1. 



If a Maass waveform vanishes in the cusp, 

lirn f(z) = 0, 

z^joo 

it is called a Maass cusp form. Maass cusp forms are square integrable over 
the fundamental domain, (/, /) < oo, where 

(/= 9) = / fg 
J r\n 

is the Petersson scalar product. 

According to the Roelcke-Selberg spectral resolution of the Laplacian jlfil 
117) . its spectrum contains both a discrete and a continuous part. The discrete 
part is spanned by the constant eigenfunction /o and a countable number of 
Maass cusp forms /i, /2, /3, ■ ■ ■ which we take to be ordered with increasing 
eigenvalues, = A < Ai < A2 < A3 < . . .. The continuous part of the 
spectrum A > 1 is spanned by the Eisenstein series E(x, 1 + ir) which are 
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known analytically |261 I27| . The Fourier coefficients of the functions Ax (1 + 
ir)E(x, 1 + ir) are given by 



b = A K (l + ir), bi = A K (l-ir), a f3 = 2 ^ |-|"\ 

A,/*£Z[i] ^ 

where 

A K {s) =An- s r{s)Q K (s) 

has an analytic continuation into the complex plane except for a pole at s = 1 . 
Normalizing the Maass cusp forms according to 

(/m fn) lj 

we can expand any square integrable function 4> G L 2 (r\TC) in terms of Maass 
waveforms, |28| . 

<t>{z) = Y J (h,<t>)Uz) + ^- f _ (E(; S ),<t>)E(z, S )ds. 
n>0 "'^ s - 1 

The eigenvalues and their associated Maass cusp forms are not known 
analytically. Thus, one has to approximate them numerically. Previous calcu- 
lations of eigenvalues for the Picard group can be found in [23 EEE EH] • 
By making use of the Hecke operators j2H EH] and the multiplicative relations 
among the coefficients, Steil |32| obtained a non-linear system of equations 
which allowed him to compute 2545 consecutive eigenvalues. We extend these 
computations with the use of HejhaPs algorithm pH], 



4 Hejhal's algorithm 

Hejhal found a linear stable algorithm for computing Maass waveforms to- 
gether with their eigenvalues which he used for groups acting on the two- 
dimensional hyperbolic plane [32], see also |3S][35]. We make use of this algo- 
rithm which is based on the Fourier expansion and the automorphy condition. 
We apply it for the Picard group acting on the three-dimensional hyperbolic 
space. For the Picard group no small eigenvalues 0<A = r 2 + l<l exist 
|37| . Therefore, r is real and the term u(y) in the Fourier expansion of Maass 
cusp forms vanishes. Due to the exponential decay of the K-Bessel function 
for large arguments (|12|l and the polynomial bound of the coefficients [3S] , 

ap = 0{\p\), Ifl-foo, 



the absolutely convergent Fourier expansion can be truncated, 
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/(*)= J2 ^yK ir (2n\(3\y)e 2 ^ +[[e\], (2) 
/3ez[i]-{o} 

\0\<M 

if we bound y from below. Given e > 0, r, and y, we determine the smallest 
M = M(e, r, y) such that the inequalities 

2-rrMy > r and K ir (2irMy) < emax(K ir (x)) 

X 

hold. Larger y allow smaller M. In all truncated terms, 



[e}}= J2 a yK ir (^\P\y)^ m0x , 



/3ez[i]-{o} 

\P\>M 

the K-Bessel function decays exponentially in |/3|, and already the K-Bessel 
function of the first truncated summand is smaller than e times most of the K- 
Bessel functions in the sum of Thus, the error [[e]] does at most marginally 
exceed e. The reason why [[e]] can exceed e somewhat is due to the possibility 
that the summands in © cancel each other, or that the coefficients in the 
truncated terms are larger than in J5Jl. By a finite two-dimensional Fourier 
transformation the Fourier expansion J5| is solved for its coefficients 

a 7 yAV(27r| 7 | y ) = — ^ £ f(x + jy) e~ 2 ^ +[[s}}, (3) 
1 ^' xexfi] 



where X[i] is a two-dimensional equally distributed set of (2Q) 2 numbers, 

X H -{ fc °^, lfcl ; *i = -Q + i-Q + l,..-,Q-|,Q-2. ^M}, 



with 2Q > M + |7|. 

By automorphy we have 

/(*) = /(**)> 

where z* is the _T-pullback of the point z into the fundamental domain J- ', 

z* = jz, 7 e r, z* e T . 

Thus, a Maass cusp form can be approximated by 

f(x+jy) = f(x*+jy*)= a0y*K ir (2ir\(3\y*)e 2 * m ^ +[{e}), (4) 

/3ezp]-{0} 

\P\<M 

where y* is always larger or equal than the height yo of the lowest points of 
the fundamental domain 
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allowing us to replace M(e, r, y) by Mo = M(e, r, yo). 

Choosing y smaller than yo the -T-pullback z i— > z* of any point into the 
fundamental domain J 7 makes at least once use of the inversion z i— » — z -1 , 
possibly together with the translations z^z+1 and z i— » z + i. This is called 
implicit automorphy, since it guarantees the invariance f(z) = /(— z -1 ). The 
conditions /(z) = f(z+l) and /(z) = f(z + i) are automatically satisfied due 
to the Fourier expansion. 

Making use of the implicit automorphy by replacing f(x + jy) in @ with 
the right-hand side of |@} gives 



a^K^-K^y) 



(20) 2 ^ 

V ^' iEX[i] /3eZ[i]-{0} 
|/3|<M 



^ 0^*^(2^/%*) e 



27ri5R/3x* -27riSR7a; 



-[[2e] 



(5) 



which is the central identity in the algorithm. 

The symmetry in the Picard group and the symmetries of the fundamental 
domain imply that the Maass waveforms fall into four symmetry classes ;3.2; 
named D, G, C, and H, satisfying 



D 

G 
C 
H 



f(x+jy) = f(ix+jy) = f(-x+jy), 
f{x + jy) = f{ix + jy) = -f(-x + jy), 
f(x + jy) = -ffa + iv) = f(-x+jy), 
f(x+jy) = -f{ix+]y) = -f(-x+]y), 



respectively, see figure |$l from which the symmetry relations among the coef- 
ficients follow, 



D 

G 
C 
H 



a/3 = <k/3 = 

a>0 = (Hp = 



l /3' 



ap 
ap 



-(hp 

-Clip 



l /3' 



Defining 



where s a p is given by 



creS/3 



SapClp 



and 
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XQ 




XQ 



XO 




y 



x 




Fig. 8. The symmetries D, G, C, and H from top left to bottom right. 



(7 G S/3 = 



{0, -0, -i0, 0, i0, -0, -i0} \i~0i {0, i0, -0, -i0}, 



^{0,10, -0,-10} else, 

the Fourier expansion Q of the Maass waveforms can be written 

f(z) = u(y)+ a p yK ir (2ir\0\y)cs{0,x), 
/3ez[i]-{o} 

where the tilde operator on a set of numbers is defined such that 
Id, (J $ K = X, and P| = 



.rex 



holds. 



Forgetting about the error [[2e]] the set of equations (JSJ) can be written as 



J2 Vy (r,y)a = O, 7 e Z[i] - {0}, | 7 | < M , 



(6) 



pez\i]-{o} 
\P\<m q 



where the matrix V = (Vyp) is given by 

V 10 {r, y) = #{-7 e S 7 }yK ir (27r| 7 |y)^ - £ y* K ir (2ir\0\y*) cs(0, x*) cs( 7> -x). 
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Since y < yo can always be chosen such that Ki r (2ir\"f\y) is not too small, the 
diagonal terms in the matrix V do not vanish for large [7] and the matrix is 
well conditioned. 

We are now looking for the non-trivial solutions of @ for 1 < I7I < Mq 
that simultaneously give the eigenvalues A = r 2 + 1 and the coefficients ap. 
Trivial solutions are avoided by setting one of the coefficients equal to one, 
a a = l. Here we choose a to be 1, 2 + i, 1, and 1 + i, for the symmetry classes 
D, G, C, and H, respectively. 

Since the eigenvalues are unknown we discretize the r axis and solve for 
each r value on this grid the inhomogeneous system of equations 

^ 7 / 3 (r,y #1 )a /3 = -T/ 7Ct (r,y# 1 ), 1 < | 7 | < M , (7) 

/3SZ[i]-{0,a} 
\0\<Mo 

where y #1 < y is chosen such that K ir (2Tr\j\y : ^ 1 ) is not too small for 1 < 
I7I < Mo- A good value to try for y^ 1 is given by 2TrM y# 1 = r. 
It is important to check whether 

9y= Vyf3(r,y* 2 )af3, 1 < M < M), 

/3ez[i]-{o} 

\/3\<M 

vanishes where y #2 is another y value independent of Only if all g 1 vanish 
simultaneously the solution of 0) is independent of y. In this case A = r 2 + 1 
is an eigenvalue and the a^'s are the coefficients of the Fourier expansion of 
the corresponding Maass cusp form. 

The probability to find an r value such that all g~ l vanish simultaneously is 
zero, because the discrete eigenvalues are of measure zero in the real numbers. 
Therefore, we make use of the intermediate value theorem where we look for 
simultaneous sign changes in g 1 . Once we have found them in at least half 
of the g 7 's we have found an interval which contains an eigenvalue with high 
probability. By some bisection and interpolation we can see if this interval 
really contains an eigenvalue, and by nesting up the interval until its size 
tends to zero we obtain the eigenvalue. 

In order not to miss eigenvalues which lie close together nor to waste CPU 
time with a too fine grid, we use the adaptive r grid introduced in |38j . 

5 Eigenvalues for the Picard group 

We have found 13950 eigenvalues of the Laplacian for the Picard group in the 
interval l<A = r 2 + l< 19601. 4115 of them belong to eigenfunctions of 
the symmetry class D, 2805 to G, 3715 to C, and 3315 to H. The smallest 
eigenvalue is A = r 2 + 1 with r — 6.6221193402528 which is in agreement 
with the lower bound A > =|- ^7\. Table shows the first few eigenvalues of 
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Table 1. The first few eigenvalues of the Laplacian for the Picard group. Listed is 
r. related to the eigenvalues via A = r 2 + 1. 



D G C H 



8.55525104 
11.10856737 
12.86991062 
14.07966049 
15.34827764 
15.89184204 
17.33640443 
17.45131992 
17.77664065 
19.06739052 
19.22290266 
19.41119126 
20.00754583 
20.70798880 
20.81526852 
21.42887079 
22.12230276 
22.63055256 
22.96230105 
23.49617692 
23.52784503 
23.88978413 
24.34601664 
24.57501426 
24.70045917 
25.47067539 
25.50724616 
25.72392169 



17.45131992 



20.70798880 



22.96230105 



23.88978413 



25.72392169 



6.62211934 
10.18079978 
12.11527484 
12.87936900 
14.14833073 
14.95244267 
16.20759420 
16.99496892 
17.86305643 
18.24391070 
18.83298996 
19.43054310 
20.30030720 
20.60686743 
21.37966055 
21.44245892 
21.83248972 
22.58475297 
22.85429195 
23.49768305 
23.84275866 
23.89515755 
24.42133829 
25.03278076 
25.42905483 
25.77588591 
26.03903968 
26.12361823 



12.11527484 

14.95244267 

16.99496892 
17.86305643 

19.43054310 
20.30030720 

21.37966055 

21.83248972 
22.58475297 

23.49768305 

23.89515755 
24.42133829 
25.03278076 

25.77588591 

26.12361823 
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each symmetry class. They agree with those of Steil [32] up to five decimal 
places. We next regard the statistics of the eigenvalues. First, we compare the 
output of our algorithm with Weyl's law and higher order corrections drawn 
from [3J||. This serves as a check whether we have found all eigenvalues. We 
then find it necessary to correct one of the terms in numerically. Finally, 
we regard the spectral fluctuations and find that the nearest-neighbor spacing 
distribution closely resembles that of a Poisson random process as predicted 
by and previously observed by [321 - 

In the first step we consider the level counting function 

AT(r) = #{i | r t <r} 

and split it into two parts 

N(r) = N(r)+N fluc (r). 

Here N is a smooth function describing the average increase in the number of 
levels and Nfi uc describes the fluctuations around the mean such that 



1 



R 



lim - / N fluc (r)dr = 0. 

R— *OC It J i 

The average increase in the number of levels is given by Weyl's law 0U]^2 an d 
higher order corrections have been calculated by Matthies [3J5]. She obtained 

N(r) = ^-r 3 + a 2 r\ogr + a 3 r + a 4 (8) 

with the constants 



fl 2 - -|?> 



= i[§log2 + f logTT - lo g r(i) + § log(2 + y/S) 



«3 



a 4 - - 5 . 



We compare our results for N(r) with (JSJ) by defining 

Nfiuc(r) = N{r) - N(r). (9) 

Nfi uc fluctuates around zero or a negative integer whose absolute value gives 
the number of missing eigenvalues, see figure El Unfortunately, our algorithm 
does not find all eigenvalues in one single run. In the first run it finds about 
97% of the eigenvalues. Apart from very few exceptions the remaining eigen- 
values are found in the third run. To be more specific, we plotted Nfi uc de- 
creased by h, because N(r) — N(r) is approximately \ whenever A = r 2 + 1 
is an eigenvalue. Furthermore, we took the eigenvalue A = into account. 
We remark that we never find more eigenvalues than predicted by JSj. A plot 
indicating that Nfi uc fluctuates around zero is shown in figure 1101 where we 
plotted the integral 
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Nfluc 




Fig. 9. Nfi uc (ri) as a function of i fluctuating around zero. 



Fig. 10. 7 as a function of R showing that / 



I(R) = i j* N fluc (r)dr. 



(10) 



So far, everything seems to be consistent. Taking the desymmetrized spectra 
into account (JSJ is modified 



N(r) 



247T 



r'+bir + b 2 r log r + b 3 r + b 4 



(11) 



with the constants depending on the symmetry class. For the symmetry class 
D the constants are given in |39| as 

h = 24, 
h — 13 

&3 = ^[-ilog2+f logTr-logr(i) 

+ § log(2 + v^) + | log(3 + 2V2) + f ], 
f. — 47 



For G 
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For C 



h 
bi 
b s 

bi 



bi 
b 2 
b 3 
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+ | log(2 + y/3) + \ log(3 + 2V2) 

25 
' 72 ' 



J_ 

96' 
1 

8?r 

— [- 

4-7T 116 



[^log2 + |log7r-lo g r(i) 
+ § log(2 + V3) - i log(3 + 2V2) 



125 
576 ' 



And for H 



96' 
1_ 

87T ' 

-L[§log2 + flog7r-lo g r(A) 
+ § log(2 + V3) - \ log(3 + 2V2) + |], 

163 
576 ' 

Let {ri} be a sequence related to the consecutive eigenvalues A = r 2 + 1. If we 
plot Nfi uc (ri) as a function of i for the desymmetrized spectra we obtain small 
deviations which can hardly be seen in figure ITT1 But if we plot the integral 
l(TT)jl we see that Nfi uc does not really fluctuate around zero. Instead, in figure 
El w e see systematic deviations, but the discrepancy is much less than one 
eigenvalue for each symmetry class. Since the number of eigenvalues is integer- 
valued we do not assume that we have found less or too many eigenvalues. 
Therefore, we fit the constants 61,62,^3,^4 in <|H[1 and obtain new constants 
for each of the symmetry classes. Since the integrals I(R) in figure IT2"1 show 
a linear behavior, the constants 61 and 62 seem to be correct. We thus only 
change the constants 63 and 64 by fitting them numerically. For the symmetry 
class D the new constants are 

63 = 0.8639... instead of 63 = 0.8679..., 

64 = -0.288... instead of 64 = -0.653... . 

For G 



h = 

62 = 

63 = 

64 = 



63 = 0.0285.. 

6 4 = -0.184. 



instead of 63 = 0.0324.. 
instead of 6 4 = —0.347. 
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2000 . 4000 

Fig. 11. Nfi uc (n) as a function of i for each symmetry class. The symmetry classes 
are D, G, C, H from top to bottom. 



50 100 p 1 50 

Fig. 12. 7 as a function of R showing the systematic deviations from I R ^°° > q 
Each curve belongs to one of the symmetry classes D, G, C, H. 

For C 

6 3 = 0.0150... instead of 63 = 0.0111..., 

b 4 = -0.062... instead of 64 = 0.217... . 

And H 

63 = 0.0702... instead of 63 = 0.0662..., 

64 = 0.034... instead of 64 = 0.283... . 

In figure El w e present the integral (|1C)|) with the corrected constants. 

Now we are able to regard the spectral fluctuations. We unfold the spec- 
trum, 



Xi = N( n ), 
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Fig. 13. J as a function of R with the corrected constants. Each curve belongs to 
one of the symmetry classes. The curves are quite indistinguishable from 0. 



in order to obtain rescaled eigenvalues x% with a unit mean density. Then 

Si — Xi-\-\ Xi 

defines the sequence of nearest-neighbor level spacings which has a mean value 
of 1 as i — * oo. We find that the spacing distribution comes close to that of a 
Poisson random process, 

P Poisson ( S ) — G , 

see figures ITU to 1 171 as opposed to that of a Gaussian orthogonal ensemble of 
random matrix theory, 

P G oe(s) ~ | se_fs2 ■ 
The integrated distribution, 

I(s) = f P(t) dt, 
Jo 

showing the fraction of spacings up to a given length is also shown in figures 
El to El The spacing distributions of the desymmetrized spectra are in ac- 
cordance with the conjecture of arithmetic quantum chaos. Also in agreement 
with the conjecture is that we have not found any degenerate eigenvalues 
within each symmetry class. But taking the eigenvalues of all four symmetry 
classes together systematic degeneracies occur due to the following: 

Theorem 1 (Steil [32|). If X = r 2 + 1 is an eigenvalue corresponding to an 
eigenfunction of the symmetry class G resp. H then there exists an eigenfunc- 
tion of the symmetry class D resp. C corresponding to the same eigenvalue. 

These degeneracies were first observed by Huntebrinker [201 and later ex- 
plained by Steil [32] with the use of the Hecke operators The Hecke 
operators are defined by 
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2 s 4 

Fig. 14. Level spacing distribution for the symmetry class D. The abscissa displays 
the spacings s. The dashed curve starting at the origin is the integrated distribution. 
For comparison, the full curves show a Poisson distribution. 




Fig. 15. Level spacing distribution for the symmetry class G. 




2 s 4 

Fig. 17. Level spacing distribution for the symmetry class H. 
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T 7ff (z) = ^ ]T g{{adr^{az + b)(d)-\ad)^), 7 S Z[i] - {0}. 

arf— 7 
f>(mod eZ) 
iftd>0, Qd>0. 

They are self adjoint operators which commute with the Laplacian and among 
each other. One can therefore simultaneously diagonalize these operators. The 
corresponding Maass cusp forms are then called Hecke eigenfunctions. The 
eigenvalue equation of the Hecke operators reads 

T 7 g r (z) = t 7 g r (z), 7 € Z[i] - {0}, 

where each Hecke eigenfunction is either identical to a Maass cusp form with 
a given symmetry or to a superposition of Maass cusp forms corresponding to 
the same eigenvalue A = r 2 + 1, but to different symmetry classes, 

9r(z) = ^ C n f n (z). 

(A+(r 2 + l))f n (z)=0 

The Hecke operators are multiplicative, 

T 1 T [j g r {z) = ^ T ^.9r(z), 

and the Hecke eigenvalues are connected to the Fourier coefficients, 

6 7 = M 7 , 7 e Z[i] - {0}, 

where the Fourier coefficients & 7 of the Hecke eigenfunctions are given by 

6 7 ^ C n (X^y n 

n 

and the index n at the Fourier coefficients of the Maass cusp forms a 7 = a 7i „ 
means that they belong to the Fourier expansion of the n-th Maass cusp form 
Uz). 

Lemma 1 (Steil [32 ). If g r (z) is a Hecke eigenfunction that does not vanish 
identically, then: 

(i) Its first Fourier coefficient is never zero, b\ ^ 0. 

(ii) A Hecke eigenfunction cannot be of symmetry class G or H. 

(Hi) Hecke eigenfunctions can always be desymmetrized such that they fall 
either into the symmetry class DUG or C U H. 

Proof (Steil's theorem). Let f n (z) be a Maass cusp form of the symmetry class 
G or H. Due to Steil's lemma, it cannot be a Hecke eigenfunction. Since one 
can diagonalize the Laplacian and the Hecke operators simultaneously, there 
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have to exist linearly independent Maass cusp forms /«.+&, k = 0,...,K 
corresponding to the same eigenvalue A = r 2 + 1 such that 

K 

^2c n+k f n+k (z) = g r (z) 

is a Hecke eigenfunction. 

At least one of these Maass cusp forms has to be of the symmetry class D 
or C in order that 

K 

b\ = Cn+k a l,n+k 
k=0 

does not vanish. 

Since the Hecke eigenfunctions can be desymmetrized such that they fall 
into the symmetry class DUG resp. C U H they are a superposition of either 
Maass cusp forms of the symmetry classes D and G or of Maass cusp forms 
of the symmetry classes C and H. Therefore, if f n {z) is of symmetry class 
G one of the / n +fe, k = 1, . . . , K is of symmetry class D, and if f n (z) is of 
symmetry class H one of the f n + k , k = 1, . . . , K is of symmetry class C. 

Based on our numerical results we now conjecture the following: 

Conjecture 4- Taking all four symmetry classes together, there are no degen- 
erate eigenvalues other than those explained by Steil's theorem. Furthermore, 
the degenerate eigenvalues which are explained by Steil's theorem occur only 
in pairs of two degenerate eigenvalues. They never occur in pairs of three or 
more degenerate eigenvalues. 

Maass cusp forms of the symmetry classes G and H indeed occur. On the 
one hand we have found a number of them numerically. On the other hand, 
Weyl's law also explains their existence. Due to this the number of eigenvalues 
whose corresponding Maass cusp forms belong to a specific symmetry class is 
in leading order independent of the choice of the symmetry class. Weyl's law 
together with Steil's theorem lead to the following: 

Conjecture 5. The sequence of non-degenerate eigenvalues in the spectrum of 
the Laplacian for the Picard group is of density zero. 

This means that as A — » oo 

#{non-degenerate eigenvalues < A} 
#{degenerate eigenvalues < A} 

Table looks as if it would contradict this conjecture. But this is due to 
the fact that only the first few eigenvalues are listed. In table |3 we list some 
consecutive large eigenvalues where we can see a better agreement with the 
conjecture. 
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Table 2. Some consecutive large eigenvalues of the Laplacian for the Picard group. 
Listed is r, related to the eigenvalues via A = r 2 + 1. 



D 




G 




c 




H 




139 


.65419675 


139 


.65419675 


139 


.66399548 


139 


,66399548 


139 


.65434417 


139 


.65434417 


139 


.66785333 


139 


.66785333 


139 


.65783548 


139 


.65783548 


139 


.66922266 


139 


.66922266 


139 


661 0.40.47 


139 


66104047 


139 


67870460 


139 


67870460 


139 


67694018 






139 


.68234200 


139 


.68234200 


139 


.68162707 


139 


.68162707 


139 


.68424704 


139 


.68424704 


139 


.68657976 






139 


.69369972 


139 


.69369972 


139 


71803029 


139 


71803029 


139 


69413379 


139 


69413379 


139 


72166907 


139 


72166906 


139 


69657741 


139. 


69657741 


139 


.78322452 


139 


.78322452 


139 


.73723373 


139 


.73723373 


139 


.81928622 


139 


.81928622 


139 


.73828541 


139 


.73828541 


139 


.81985670 


139 


.81985670 


139 


.74467774 


139 


.74467774 


139 


.82826034 


139 


.82826034 


139 


75178180 


139. 


75178180 


139 


.84250751 






139 


.75260292 


139 


.75260292 


139 


.87781072 


139 


.87781072 


139 


.79620628 


139 


.79620628 


139 


.87805540 






139 


.80138072 


139 


.80138072 


139 


88211647 


139 


88211647 


139 


.81243991 


139 


.81243991 


i on 
139 


m Toonno 
.91 (S2UU3 


i on 
139 


ni Toonno 
.91 i 82003 


i on 
139 


o i o i onoo 
.51312982 


i on 
139 


o 1 o 1 onoo 
.81312982 


139 


.91893517 






139 


.82871870 


139 


.82871870 


139 


.92397167 


139 


.92397167 


139 


.86401372 


139 


.86401372 


139 


.92721861 


139 


.92721861 


139 


.86461581 


139 


.86461581 


139 


.93117207 


139 


.93117207 


139 


.89407865 


139 


.89407865 


139 


.93149277 


139 


.93149277 


139 


.89914777 


139 


.89914777 


139 


.94067283 






139 


.90090849 


139 


.90090849 


139 


.94396890 


139 


.94396890 


139 


.91635302 


139 


.91635302 


139 


.95074070 






139 


.94071729 


139 


.94071729 


139 


.95124805 


139 


.95124805 


139 


.95080198 


139 


.95080198 


139 


.99098324 


139 


.99098324 


139 


.97043676 


139 


.97043676 
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6 Summary 

Our principal goal was to test the conjecture of arithmetic quantum chaos 
numerically in one example. For this purpose we have chosen a point particle 
moving freely in the three-dimensional and negatively curved quotient space 
of the Picard group. Identifying the solutions of the stationary Schrodinger 
equation with Maass waveforms allowed us to use Hejhal's algorithm to com- 
pute the eigenfunctions and eigenvalues numerically. Having computed 13950 
eigenvalues (and eigenfunctions), which exceeds all previous computations in 
non-integrable three-dimensional systems, we demonstrated that our numeri- 
cal results are in accordance with the conjecture of arithmetic quantum chaos. 
Within each symmetry class we do not find any degenerate eigenvalues, but 
taking all four symmetry classes together, almost all eigenvalues become de- 
generate in the limit of large eigenvalues A — > oo. This behaviour was explained 
by the interplay of the symmetries with the Hecke-operators. 
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A The K-Bessel function 

The K-Bessel function is defined by 

/■oo 

K ir (x) = / e ~ xcosht cos(ri) dt, dix > 0, r e C, 
Jo 

see Watson and is real for real arguments x and real or imaginary order 
ir. It solves the modified Bessel differential equation 

x 2 u"{x) + xu'(x) ~ (x 2 ~ r 2 )u{x) = 0, 

and decays exponentially for large arguments 

Ki r (x) ~ \j^~ e ~ X ^ or x ~~ * 00 ' (^) 

A second linearly independent solution of the modified Bessel differential equa- 
tion is the I-Bessel function 
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oo ,x\2k 



2' j-^ k\r{ir + k+ 1) 
which grows exponentially for large arguments 



- — e x for x — > oo . 
2ttx 

The K-Bessel function decreases exponentially when r increases. This can be 
compensated by multiplication with the factor e~ . 

In order to compute the K-Bessel function numerically for small or moder- 
ate imaginary order we use its continued fraction representation which follows 
from the Miller algorithm 02] ■ 
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